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We show how a general formulation of the Fluctuation-Response Relation is able to describe 
in detail the connection between response properties to external perturbations and spontaneous 
fluctuations in systems with fast and slow variables. The method is tested by using the 360- variable 
Lorenz-96 model, where slow and fast variables are coupled to one another with reciprocal feedback, 
and a simplified low dimensional system. In the Fluctuation-Response context, the influence of the 
fast dynamics on the slow dynamics relies in a non trivial behavior of a suitable quadratic response 
function. This has important consequences for the modeling of the slow dynamics in terms of a 
Langevin equation: beyond a certain intrinsic time interval even the optimal model can give just 
statistical prediction. 



I. INTRODUCTION 



One important aspect of climate dynamics is the study of the response to perturbations of the external forcings, or 
of the control parameters. In very general terms, let us consider the symbolic evolution equation: 

^-Q(X) (i.i) 

where X is the state vector for the system, and Q(X) represents complicated dynamical processes. As far as climate 
modeling is concerned, one of the most interesting properties to study is the so-called Fluctuation-Response relation 



(FRR), i.e. the possibility, at least in principle, to understand the behavior of the system (1.1) under perturbations 
(e.g. a volcanic eruption, or a change of the C O2 concentration) in terms of the knowledge obtained from its past 
time history (Leith 1975, 1978, Dymnikov and Gritsoun, 2001). 



The average effect on the variable Xi{t) of an infinitesimal perturbation S{{t) in (1.1 1, i.e. Q(X) Q(X) + Si(t), 
can be written in terms of the response matrix Rij{t). If di{t) = for < < one has: 

SXdij = V /* R^J{t - t')df^{t')dt' (1.2) 

where Rij{t) is the average response of the variable Xi at time t with respect to a perturbation of Xj at time 0. 

The basic point is, of course, how to express Rij(t) in terms of correlation functions of the unperturbed system. The 
answer to this problem is the issue of the Fluctuation-Response theory. This field has been initially developed in the 
context of equilibrium statistical mechanics of Hamiltonian systems; this generated some confusion and misleading 
ideas on its validity. As a matter of fact, it is possible to show that a generalized FRR holds under rather general 
hypotheses (Deker and Haake, 1975, Falcioni et al., 1990): the FRR is also valid in non Hamiltonian systems. It 



2 

is interesting to note that, although stochastic and deterministic systems, from a conceptual (and technical) point 
of view, are somehow rather different, the same FRR holds in both cases, see Appendix A. For this reason, in the 
following, we will not separate the two cases. In addition, a FRR holds also for not infinitesimal perturbation (Boffetta 
et al., 2003). From many aspects, the FRR issues in climate systems are rather similar to those in fluids dynamics: we 
have to deal with non Hamiltonian and non linear systems whose invariant measure is non Gaussian (Kraichnan, 2000). 
On the other hand, it is obviously impossible to model climate dynamics with equations obtained from flrst principles, 
so typically it is necessary to work with simple raw models or just to deal with experimental signals (Ditlevsen, 1999, 
Marwan et al., 2003). In addition, in climate problems (and more in general in Geophysics) the study of inflnitesimal 
perturbation is rather academic, while a much more interesting question is the relaxation of large perturbations in the 
system due to fast changes of the parameters. Numerical simulations show that, in systems with one single time scale 
(e.g. low dimensional chaotic model as the Lorenz one), the amplitude of the perturbation is not so important, (see 
Appendix A, and Boffetta et al., 2003). On the contrary, in the case of different characteristic times, the amplitude of 
the perturbation can play a major role in determining the response, because different amplitudes may affect features 
with different time properties (Boffetta et al., 2003). Starting from the seminal works of Lcith (1975, 1978), who 
proposed the use of FRR for the response of the climatic system to changes in the external forcing, many authors 
tried to apply this relation to different geophysical problems, ranging from simplified models (Bell, 1980), to general 
circulation models (North et al., 1999, Cionni et al., 2004) and to the covariance of satellite radiance spectra (Haskins 
et al., 1999). For recent works on the application of the FRR to the sensitivity problem and the predictability see 
Gritsoun and Dymnikov (1999), Gritsoun (2001), Gritsoun et al. (2002), Dymnikov and Gritsoun (2005), Dymnikov 
(2004), Abramov and Majda (2007), and Gritsoun and Branstator (2007). In most works, the FRR has been invoked 
in its Gaussian version, see below, which has been used as a kind of approximation, often without a precise idea of its 
limits of applicability. In principle, according to Lorenz (1996), one has to consider two kinds of sensitivity: to the 
initial conditions (flrst kind) and to the parameters (e.g. external forcing) of the system (second kind). On the other 
hand, if one considers just inflnitesimal perturbations, it is possible to describe the second kind problem in terms of 
the flrst one. Unfortunately, this is not true for non inflnitesimal perturbations. 

In this paper we study, in the FRR framework, systems with more than one characteristic time. In section 2 we 
recall the theoretical basis of the FRR issue. In section 3 we describe the analysis we have performed on two dynamical 
systems. The flrst one, is a model introduced by Lorenz (1996), which contains some of the relevant features of climate 
systems, i.e. the presence of fast and slow variables (see Fraedrich, 2003, for a discussion about short and long-term 
properties of complex multiscalc systems like the atmosphere). We consider, at this regard, the problem of the 
parameterization of the fast variables via a suitable renormalization of the parameters appearing in the slow dynamics 
equations, and the addition of a random forcing. The second one is a very simplified system consisting, basically, of 
a slow variable which fluctuates around two states, coupled to fast chaotic variables. The specific structure of this 
system suggests a modeling of the slow variable in terms of a stochastic differential equation. We will see how, even 
in absence of a Gaussian statistics, the correlation functions of the slow (fast) variables have, at least, a qualitative 
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resemblance with response functions to perturbations on the slow (fast) degrees of freedom. In addition, although 
the average response of a slow variable to perturbations of the fast components is zero, the influence of the fast 
dynamics on the slow dynamics cannot be neglected. This fact is well highlighted by a non trivial behavior of a 
suitable quadratic response function (Hohenberg and Shraiman, 1989). In the framework of the complexity in random 
dynamical systems, one has to deal with a similar behavior: the relevant 'complexity' of the system is obtained by 
considering the divergence of nearby trajectories evolving with two different noise realizations (Paladin et al., 1995). 
This has important consequences for the modeling of the slow dynamics in terms of a Langevin equation: beyond a 
certain intrinsic time interval (determined by the shape of the quadratic response function) even the optimal model 
can give just statistical predictions (for general discussion about the skills and the limits of predictability of climatic 
models see Cane, 2003). The conclusions and the discussion of the results obtained in this work are contained in 
section 4, while the Appendices are devoted to some technical aspects. 



For the sake of completeness we summarize here some basic results regarding the FRR (see Appendix A for technical 
details). Let us consider a dynamical system X(0) X(t) = C/*X(0) whose time evolution can even be not completely 
deterministic (e.g. stochastic differential equations), with states X belonging to a iV-dimensional vector space. We 
assume: a) the existence of an invariant probability distribution p(X), for which some "absolute continuity" type 
conditions are required (see Appendix A); b) the mixing[56j character of the system (from which its crgodicity 
follows) . 

At time t = we introduce a perturbation (5X(0) on the variable X(0). For the quantity SXi{t), in the case of an 
infinitesimal perturbation (5X(0) = {6Xi{0) ■ ■ ■ 6X^(0)) one obtains: 



In the following (()) indicates the average on the unperturbed system, while () indicates the mean value of perturbed 
quantities. The operating definition of Rij{t) in numerical simulations is the following. We perturbe the variable Xj 
at time t = to with a small perturbation of amplitude SXj(0) and then evaluate the separation component SXi(t\tQ) 
between the two trajectories X(t) and X'(t) which are integrated up to a prescribed time ti = to + At. At time t — ti, 
the variable Xj of the reference trajectory is again perturbed with the same dXj{0), and a new sample SJi.{t\ti) is 
computed and so forth. The procedure is repeated M 3> 1 times and the mean response is then given by: 



II. THEORETICAL BACKGROUND ON FRR 




(III) 



where the linear response functions (according to FRR) are 




(II.2) 




M 



Usually, in non Hamiltonian systems, the shape of p(X) is not known, therefore relation (II. 2 1 does not give a very 
detailed information. On the other hand the above relation shows that, anyway, there exists a connection between 
the mean response function R^j and some suitable correlation function, computed in the unperturbed systems. 

In the case of multivariate Gaussian distribution. In /3(X) = — | j ctij^i^j + const, where {ctij} is a positive 
symmetric matrix, the elements of the linear response matrix can be written in terms of the usual correlation functions, 
= {X,{t)Xk{0))/{X,Xk), as: 

R,j{t) = J2(^,k{x,{t)Xk{0)) . (II.3) 

One important nontrivial class of systems with a Gaussian invariant measure is the inviscid hydrodynamics [57 , where 
the Liouville theorem holds, and a quadratic invariant exists (Kraichnan, 1959, Kraichnan and Montgomery, 1980, 
Bohr et al., 1998). Sometimes in the applications, in absence of detailed information about the shape of p, formula 



(II. 3 1 is assumed to hold to some extent. Numerical studies of simplified models which mimic the chaotic behavior 



of turbulent fluids show that, since that stationary probability distribution is not Gaussian, Eq. (II. 3) does not 
hold. On the other hand, the correlation functions and the response functions have similar quantitative behavior. In 
particular, in fully developed turbulence, as one can expect on intuitive ground, one has that the times characterizing 
the responses approximate the characteristic correlation times (Biferale et al., 2002, Boffetta et al., 2003). This is 
in agreement with numerical investigation (Kraichnan, 1966) at moderate Reynolds number of the Direct Interaction 
Approximation equations, showing that, although Rii{t) is not exactly proportional to the autocorrelation function 
Cii{t), if one compares the correlation times Tc{ki) (e.g. the time after which the correlation function becomes lower 
than 1/2) and the response time Tfi(ki) (e.g. the time after which the response function becomes lower than 1/2), the 
ratio Tc{ka) /Tii{ka) remains constant through the inertial range. In the turbulence context, Xi indicates the Fourier 
component of the velocity field corresponding to a wave vector fc^. 

We would like briefly to remark a subtle point. From a rather general argument (see Appendix B), one has that all 
the (typical) correlation functions, at large time delay, have to relax to zero with the same characteristic time, related 
to spectral properties of the operator L which rules the time evolution of the probability density function P(X,i): 

-P(X,i)-LP(X,t). (II.4) 

Using this result in a blind way, one has the apparently paradoxical conclusion that, in any kind of systems, all the 
correlation functions, relative to degrees of freedom at different scales, relax to zero with the same characteristic time. 
On the contrary, in systems with many different characteristic times (e.g. fully developed turbulence), one expects a 
whole hierarchy of times distinguishing the behavior at different scales (Frisch, 1995). The paradox is, of course, only 
apparent since the above argument is valid just at very long times, i.e. much longer than the longest characteristic 
time, and therefore, in systems with fluctuations over many different time-scales, this is not very helpful. 
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III. RESPONSE OF FAST AND SLOW VARIABLES 

Systems with a large number of components and/or with many time scales, e.g. climate dynamics models, present 
clear practical difficulties if one wants to understand their behavior in detail. Even using modern supercomputers, it 
is not possible to simulate all the relevant scales of the climate dynamics, which involves processes with characteristic 
times ranging from days (atmosphere) to 10^-10'^ years (deep ocean and ice shields), see Majda et al. (2005) and 
Majda and Wang (2006). 

For the sake of simplicity, we consider the case in which the state variables evolve over two very different time 
scales: 



dt 

dJ^^l 
dt e 



f(X„X/) (III.l) 
g(X.,X/) (III.2) 



where Xs and X/ indicate the slow and fast state vectors, respectively, e ^ 1 is the ratio between fast and slow 
characteristic times, and both f and g are 0(1). A rather general issue is to understand the role of the fast variables 
in the slow dynamics. From the practical point of view, one basic question is to derive effective equations for the slow 
variables, e.g. the climatic observable, in which the effects of the fast variables, e.g. high frequency forcings, are taken 
into account by means of stochastic parameterization. Under rather general conditions (Givon et al., 2004), one has 
the result that, in the limit of small e, the slow dynamics is ruled by a Langevin equation with multiplicative noise: 

dX., 



dt 



= fe//(X,) + a(X3)77 (III.3) 



where i] is a white-noise vector, i.e. its components are Gaussian processes such that {rjiit)) and {'qi{t)rjj(t')) ~ 
Sij6{t — t'). Although there exist general mathematical results (Givon et al., 2004) on the possibility to derive eq. 



(III. 3 1 from (III.l ) and (III. 2 1, in practice one has to invoke (rather crude) approximations based on physical intuition 
to determine the shape of f^jj and a (Mazzino et al., 2005). At this regard, see also the contribution to the volume 
by Imkeller and von Storch (2001) about stochastic climate models. For a more rigorous approach in some climate 
problems see Majda et al. (1999, 2001) and Majda and Franzke (2006). 

In the following, we analyse and discuss two models which, in spite of their apparent simplicity, contain the basic 
features, and the same difficulties, of the general multiscale approach: the Lorenz-96 model (Lorenz 1996) and a 
double-well potential with deterministic chaotic forcing. 

A. The Lorenz-96 model 



First, let us consider the Lorenz-96 system (Lorenz 1996), introduced as a simplified model for the atmospheric 
circulation. Define the set {xk{t)}, for k — l,...,Nk, and {yk,j{t)}, for j — l,...,Nj, as the slow large-scale variables 



and the fast small-scale variables, respectively (being Nk = 36 and Nj = 10). Roughly speaking, the {x^j's represent 
the synoptic scales while the {y/cj l's represent the convective scales. The forced dissipative equations of motion are: 

AT,- 

~ — T.u 1 (t.i. n — T.l. Ill — IJT.r. A- V -\- r-i 

dt 



'^^'^ --Xk-i{xk-2- Xk+i) - vxk + F + ci'^Uk^j (III. 4) 

= -cbykj+i{yk,j+2 - Vkj-i) - ciyyk,j + ciXk (III. 5) 

where: = 10 is the forcing term, = 1 is the linear damping coefficient, c = 10 is the ratio between slow and fast 
characteristic times, 6 = 10 is the relative amplitude between large scale and small scale variables, and ci = c/b = 1 
is the coupling constant that determines the amount of reciprocal feedback. 

Let us consider, ffi-st, the response properties of fast and slow variables, see Figs. [l]and[2] 

In Fig. [l] the autocorrelation Cjj{t) and self-response Rjj{t) refer to the fast variable yk.j{t), with fixed k and j. 
It is well evident how, even in absence of a precise agreement between autocorrelations and self-response functions 
(due to the non Gaussian character of the system), one has that the correlation of the slow (fast) variables have at 
least a qualitative resemblance with the response of the slow (fast) variables themselves. 

The structure of the Lorenz-96 model includes a rather natural set of quantities that suggests how to parameterize 
the effects of the fast variables on the slow variables, for each k. Let us indicate with — X^j^i Vk.j the term 
containing all the Nj fast terms in the equations for the Nk slow modes. In the following, we will see that, replacing 
the deterministic terms {zfej's in the equations for the {x^j's with suitable stochastic processes, one obtains an 
effective model able to reproduce the main statistical features of the slow components of the original system. 

It's worth- noting, from Fig. |3] that Ckk{t) and Cz^{t), the autocorrelation of the cumulative variable Zk{t), are 
rather close to each other. This suggests that Zk{t) must be correlated to Xk{t), in other words, the cumulative effects 
of the Nj fast variables yk.j{t) on Xk.(t) are equivalent to an effective slow term, proportional to Xk{t). 

We look, therefore, for a conditional white noise parameterization that takes into account this important information 
given by the structure of the Lorenz-96 model equations. Let us write the effective equations for the slow modes as 

— = -Xk-i{xk~2 - Xk+i) - {v + v')xk + {F + F') + C2-r]k (III. 6) 

where 77^, are uncorrelated and normalized white-noise terms. Some authors, Majda et al (1999, 2001) and Majda 
and Franzke (2006), using multiscale methods, have obtained effective Langevin equations for the slow variables of 
systems having the same structure as the Lorenz-96 model. 

Basically we can say that, in the effective model for the slow variables, one parameterizes the effects of the fast 
variables with a suitable renormalization of the forcing, F F + F', of the viscosity, u ^ 1^+ i'' , and the addition of 



a random term. In other words, we replace the Zk — X^^i Vk.j terms in (III. 4) with stochastic processes z^ depending 



on the slow variables Xk'- 

dxk 



-Xk-iixk-2 - Xk+i) - vxk + F + ciZk (III. 7) 



where 



Zk = — {-v'xk + F' + C2?7fe) 

Cl 



(III, 



with Cl is a new couphng constant. We notice that eq. (III. 7 1 has the same form of eq. (III.4I. With a proper choice 



of u' , F' and C2 in (III. 6), v' = —0.3, F' ~ 0.25, C2 = 0.3, which imphes ci = 0.25 in ( |III.7[ ), one can reproduce 
the statistics of Xk and Zk to a very good extent, see at this regard Fig. [4] Fig. [5] and Fig. [6] Of course the above 
described parameterization of the fast variables is inspired to the general 'philosophy' of the Large-Eddy Simulation 
of turbulent geophysical flows at high Reynolds numbers (Moeng, 1984, Moeng and Sullivan, 1994, Sullivan et al., 
1994). 

The FR properties of the stochastic Lorenz-96 slow variables are reported in Fig. [Tj 

Let us come back to the response problem. Of course the mean response of a slow variable to a perturbation on 
a fast variable is zero. However, this does not mean that the effect of the fast variables on the slow dynamics is not 
statistically relevant. Let us introduce the quadratic response of Xk{t) with respect to an infinitesimal perturbation 
on ykj{0), for fixed k and j: 

1/2 



Sxkitf 



^J/fc,i(0) 



(III.9) 



Considered that in all simulations the initial impulsive perturbations on the ykj is kept constant, 5yk.j{0) ~ A, with 



A ^ iUk it is convenient to take the average of (III. 9 1 over all j's, at a fixed k, and introduce the quantity 



(III. 10) 



where with s and / we label the slow and fast variables, respectively. In the case of the Lorenz-96 system, all the ykj 
variables, at fixed fc, are statistically equivalent, and have identical coupling with Xk, so that R^!'^ (t)/ A coincides with 
^fcf We report in Fig. jsjthe behavior of R'f^{t), for both ( |lll.4[ ) and |lll.7| . As regards to the stochastic model, 
the analogous of ( Ill.lOp is defined as follows. One studies the evolution of Sxk{t) as difference of two trajectories 
obtained with two different realizations of the {ry^j's. It is worth stressing that the behavior of 5xk(t) under two noise 
realizations can be very different from the behavior of Sxk(t) under the same noise realization (see Appendix C). This 
aspect will be considered again in the next section. 



B. A simplified model 



In order to grasp the essence of systems with fast and slow variables, we discuss now a toy climate model in which 
the 'climatic' variable fluctuates between two states. Consider a four dimensional state vector q = (90,91,92,93) 
whose evolution is given by: 

dqo 



dt 



cqi 



(III.ll) 



dqi 1. , 
at e 



dq2 
dt 



;[-9i93 + ''L^i - 92] 



If , 1 



(III. 12) 
(III. 13) 
(III. 14) 



This four equation system will be named the deterministic DW model. The subsystem formed by (III. 12 1, (III. 13 1 



and (III. 14 1 is nothing but the well-known Lorenz-63 model (Lorenz 1963), in which the constant ?has the function 
of rescaling the characteristic time. In absence of coupling (c = 0) between go said qi, the unforced motion equation 
holds for the slow variable x — qg: 



dx 
~di 



dV_ 
dx 



2^/Hx - x^ with V{x) = H - \/Hx 



1 



(III. 15) 



The system (III. 15 1 has one unstable steady state in a;o = corresponding to the top of the hill of height H, and two 
stable steady states in xi/2 — ±(4i/)^/*, i.e. the bottom of the valleys. The presence of the coupling (c ^ 0) between 



slow and fast variables can induce transitions between the two valleys. The parameters in (III. 11), (III. 12 1, (III. 13 1 



and (III. 14) are fixed to the following values: ctl ~ 10, = 28, &l = 8/3, i.e. the classical set-up corresponding to 
the chaotic regime for the Lorenz-63 system; H ^ 4, the height of the barrier; c — 0.5, the coupling constant that 
rules the transition time scale of <7o(i) between the two valleys; by setting e = 1, the ratio e between fast and slow 



characteristic times, see (III.ll and (III.2I, is O(10~^). 

Since the time scale of the qo{t) well-to-well transitions may be considerably longer, depending on 1/c, than the 
characteristic time of qi{t), of order 0(1), we refer to go ^ the slow variable, or the low-frequency observable, and 
to qi as the fast variable, or the high-frequency forcing, of the deterministic DW model. It can be easily shown that, 
for c — 0, small perturbations Ago around the two potential minima at ±(4iJ)^/"' relax exponentially to zero with 
characteristic time l/'iy/H. For sufficiently large values of c, the climatic variable qQ{t) jumps aperiodically back and 
forth between the two valleys, driven by the chaotic signal qi{t), see Fig. [9] and Fig. [lO] 



The main statistical quantities investigated to analyse the DW model are the following: 

a) the probability density function of the slow variable qo', 

b) the probability density function of the well-to-well transition time t^, p{te)\ 

c) the slow and fast auto-correlation functions (ACF) Cii{t) = {qi{t)qi{0)) / (qf), with i = 0, 1; 



d) the slow and fast self-response functions (ARF) Ru^t) = 6qi(t)/Sqi(0), with i = 0,1; 

e) the quadratic cross-response function of the slow variable qo{t) with respect to the fast variable gi(0). 

Of course i?oi(t), i.e. the mean response of qoit) to a perturbation on qi{0), is zero for trivial symmetry arguments. 
On the other hand, the quadratic response: 



Sqoitr 


1/2 


SqiiO 





(III.16) 
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can give relevant physical information. Even in this case, since in all simulations the initial perturbation on qi{0) is 
kept constant, Sqi{0) = A <C {qf)^^^, it is convenient to define as mean quadratic response of the slow variable (s) 
with respect to the fast variable (/) the quantity = A • R'^^it). The long-time saturation level of R'fj!(t) is of 

the order of the distance between the two climatic states. 

With the current set-up, slow and fast variable have characteristic times which differ by an order of magnitude 
from each other, while the statistics of go is strongly non Gaussian. Because of the skew structure of the system, 
i.e. the fast dynamics drives the slow dynamics but without counter-feedback, one expects that, at the least in the 
limit of large time scale separation, the joint PDF can be factorized, with an asymptotic PDF for go of the form 
po ^ K ■ e^^"^-'' where K is a, normalization constant. 

The FR properties of the deterministic DW model, for the fast and slow variables, are shown in Fig. [TT] and Fig. 
[T2| respectively. 

The slow self-response i?oo(i) initially decreases exponentially with characteristic time {H = 4), i.e. the 

same behavior of the relaxation of a small perturbation near the bottom of a valley for c = 0. Then, i?oo(^) relaxes 
to zero much more slowly. It is natural to assume that this is due to the long-time jumps between the valleys. It 
is well evident that Rqq behaves rather differently from Coo, while Rn and Cn have, at least, the same qualitative 
shape. On the other hand, the autocorrelation (self-response) time scales of the two variables differ from each other 
of a factor ^10, compatibly with the fact that the ratio between fast and slow characteristic times is e ~ 0.1, for the 
current set-up (e = 1). 

Since the statistics is far from being Gaussian, the 'correct' correlation function which satisfies the FR theorem, for 
the slow variable, has the form: 

dpe{qo,qi,q2,q3) 



C{t) = - qoit) 



t=0 



(III. 17) 



dqa 

where Pe{qa, qi,q2, qs) is the (unknown) joint PDF of the state variable of the system at a fixed e. In the limit of large 
time separation, i.e. for e^— > 0, one expects that the asymptotic PDF /3o('?Oi 9ij 92, 93) is factorized: 

^0(90,91,92,93) =ife-^'="(*VL('?i,'72,g3) (111.18) 

where A" is a normalization constant, and is the PDF of the Lorenz-63 state variable. Under this condition, the 
right correlation function predicted by the FRR has a relatively simple form: 

dVeffiqa) 



C{t) = qoit) 



(III. 19) 



^90 

where Vg// indicates the effective potential. For e ^ 10~^ (corresponding to = 1) we have checked numerically 
that the joint PDF is not yet factorized, while for a very small ratio between the characteristic times, e ~ 10~^ 
(corresponding to e= 10^^), the form (III. 18 1 holds and, taking V^ff (x V, we obtain a very good agreement between 



Rooit) and C(t), see Fig. 13 



The cross-response properties of the DW model, measured by the quantity r[,''} (t), are reported in Fig. 



18 



We will 



consider again later this issue when discussing the stochastic modeling. While the mean (slow-to-fast) cross-response 
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is null (not shown), its fluctuations grow with time. This means that an initial uncertainty on the fast variables has 
consequences for the predictability of the slow variable, since it induces a mean separation growth between two initially 
close 'climatic' states of the go variable. At small times, grows exponentially in time, i.e. it is driven by the 

chaotic character of the fast variable while, at very long times, the well-to-well aperiodic jumps play the dominant 
role and the growth speed eventually decreases to zero until saturation sets in. 

Let us now consider a stochastic model for the slow variable qoit), obtained by replacing the fast variable qi, in the 
equation for qq, with a white noise. One has a Langevin equation of the kind: 



(t) = 2VHqoit) - go' + ^ ■ m (111.20) 



where ^{t) is a Gaussian process with {^{t)) = and {^{t)^{t')) = S{t - t'). We call eq. (III.20I the WNDW model 



The value a — 19.75 is determined by requiring that the PDFs of the well-to-well transition times have the same 



asymptotic behavior (i.e. exponential tail with the same exponent), see Fig. 14 

Let us notice that, in this case, because of the skew structure of the original system, the stochastic modeling is 
(relatively) simple and, differently from the generic case, the noise is additive. The time signal go(^) obtained from 
the WNWD model is reported in Fig. [Tsj One observes strong similarities in the long-time transition statistics with 
respect to the deterministic model, even though the PDFs of the slow variable are quite different from one another, 
see Fig. [16] 



The FR properties of the WNDW model are reported in Fig. 17 The slow variable is distributed according to 
^-v(qo)/K ^ with K = and the FR theorem prediction is verified, i.e. one has a good agreement between 

Roo{t) and the correlation function C{t). 

We redefine, as already seen when discussing the stochastic model approximating the Lorenz-96 system, the 



quadratic cross-response function r[,''} (t) as the root mean square growth of the error Sqo{t) induced by two dif 



ferent noise realizations. 



In Fig. 18 the behavior of R^^'^) {t) for the deterministic DW system and its stochastic model is reported. The WNDW 



model is not able to reproduce the two-time behavior of the deterministic model, mainly due to the impossibility to 
control the amplitude of the initial perturbation. Because of that, the error on the climatic state of the system 
saturate very quickly, as soon as the trajectory starts jumping between the wells. 



IV. DISCUSSION AND CONCLUSIVE REMARKS 



In this paper we have presented a detailed investigation of the Fluctuation-Response properties of chaotic systems 
with fast and slow dynamics. The numerical study has been performed on two models, namely the 360-variable 
Lorenz-96 system, with reciprocal feedback between fast and slow variables, and a simplified low dimensional system, 
both of which are able to capture the main features, and related difficulties, typical of the multiscale systems. The first 
point we wish to emphasize is how, even in non Hamiltonian systems, a generalized Fluctuation-Response Relation 
(FRR) holds. This allows for a link between the average relaxation of perturbations and the statistical properties 
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(correlation functions) of the unperturbed system. Although one has non Gaussian statistics, the correlation functions 
of the slow (fast) variables have at least a qualitative resemblance with the response functions to perturbations on the 
slow (fast) degrees of freedom. The average response function of a slow variable to perturbations of the fast degrees of 
freedom is zero, nevertheless the impact of the fast dynamics on the slowly varying components cannot be neglected. 
This fact is clearly highlighted by the behavior of a suitable quadratic response function. Such a phenomenon, which 
can be regarded as a sort of sensitivity of the slow variables to variations of the fast components, has an important 
consequence for the modeling of the slow dynamics in terms of a Langevin equation. Even an optimal model (i.e. 
able to mimic autocorrelation and self- response of the slow variable), beyond a certain intrinsic time interval, can 
give just statistical predictions, in the sense that, at most, one can hope to have an agreement among the statistical 
features of system and model. In stochastic dynamical systems, one has to deal with a similar behavior: the relevant 
'complexity' of the systems is obtained by considering the divergence of nearby trajectories evolving under two different 
noise realizations. Therefore a good model for the slow dynamics (e.g. a Langevin equation) must show a sensitivity 
to the noise. 



In this Appendix we give a derivation, under general rather hypothesis, of a generalized FRR. Consider a dynamical 
system x(0) x(i) = [/*x(0) with states x belonging to a A^-dimensional vector space. For the sake of generality, we 
will consider the case in which the time evolution can also be not completely deterministic {e.g. stochastic differential 
equations). We assume the existence of an invariant probability distribution p(x), for which some "absolute continuity" 
type conditions are required (see later), and the mixing character of the system (from which its ergodicity follows). 
Note that no assumption is made on N. 

Our aim is to express the average response of a generic observable ^ to a perturbation, in terms of suitable correlation 
functions, computed according to the invariant measure of the unperturbed system. At the first step we study the 
behavior of one component of x, say Xi, when the system, described by p(x), is subjected to an initial (non-random) 
perturbation such that x(0) — > x(0) + Axg. This instantaneous kick[58 modifies the density of the system into /9'(x), 
related to the invariant distribution by p'(x) = /3(x — Axo). We introduce the probability of transition from xo at time 
to X at time t, VF(xo,0 x,t). For a deterministic system, with evolution law x(t) = L/*x(0), the probability of 
transition reduces to W{:x.o, ^ x, t) = (5(x — ?7*xo), where S{-) is the Dirac's delta. Then we can write an expression 
for the mean value of the variable Xi, computed with the density of the perturbed system: 



V. APPENDIX A: GENERALIZED FRR 




(V.l) 



The mean value of Xi during the unperturbed evolution can be written in a similar way: 




(V.2) 
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Therefore, defining Sxi 



(xi), we have: 



6xt{t) = y y i^(xo, Axo) p(xo)VF(xo,0 — > x,t)dxrfxo 
[x,{t) F(xo,Axo)^ 



where 



F(xo,Axo) 



p(xo - Axq) - p(xo) 
p(xo) 



(V.3) 



(V.4) 



Let us note liere that the mixing property of the system is required so that the decay to zero of the time-correlation 
functions assures the switching off of tire deviations from equilibrium. 

For an infinitesimal perturbation (5x(0) = {Sxi{0) ■ ■ ■ Sxn{0)), if p(x) is non- vanishing and differentiable, the function 



in (V.4 1 can be expanded to first order and one obtains: 



din p(x.) 



dxj 



>fa,(0) 



t=o/ 



Y^R,j{t)5x,{Q) 



(V.5) 



which defines the linear response 



R,j{t)^~{x,{t) 



din p{x) 
dxi 



(V.6) 



of the variable Xi with respect to a perturbation of Xj. One can easily repeat the computation for a generic observable 
A(x): 



SA{t) 



)6x,{0) 



(V.7) 



For Langevin equations, the differentiability of p(X) is well established. On the contrary, one could argue that in 
a chaotic deterministic dissipative system the above machinery cannot be applied, because the invariant measure is 
not smooth at all. Typically the invariant measure of a chaotic attractor has a multifractal character and its Renyi 
dimensions dq are not constant (Paladin and Vulpiani, 1987). In chaotic dissipative systems the invariant measure is 
singular, however the previous derivation of the FRR is still valid if one considers perturbations along the expanding 
directions. For a mathematically oriented presentation see Ruelle (1998). A general response function has two 
contributions, corresponding respectively to the expanding (unstable) and the contracting (stable) directions of the 
dynamics. The first contribution can be associated to some correlation function of the dynamics on the attractor (i.e. 
the unperturbed system) . On the contrary this is not true for the second contribution (from the contracting directions) , 
this part to the response is very difficult to extract numerically (Cessac and Sepulchre, 2007). In chaotic deterministic 
systems, in order to have a differentiable invariant measure, one has to invoke the stochastic regularization (Zeeman 
1990). If such a method is not feasible, one can use the direct approach by Abramov and Majda (2007). For a study 
of the FRR in chaotic atmospheric systems, see Dymnikov and Gritsoun (2005) and Gritsoun and Branstator (2007). 
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Let us notice that a small amount of noise, that is always present in a physical system, smoothen the p{'x.) and the 
FRR can be derived. We recall that this "beneficial" noise has the important role of selecting the natural measure, and, 
in the numerical experiments, it is provided by the round-off errors of the computer. We stress that the assumption 
on the smoothness of the invariant measure allows to avoid subtle technical difficulties. 



VI. APPENDIX B: A GENERAL REMARK ON THE DECAY OF CORRELATION FUNCTIONS 



Using some general arguments one has that all the (typical) correlation functions at large time delay have to relax to 
zero with the same characteristic time, related to spectral properties of the operator L which rules the time evolution 
of the P(X,i): 

^P(X,^)=LP(X,<). (VI.l) 

In the case of ordinary differential equations 

dXi/dt= Qi(X.) i^l,---,N (VI.2) 

the operator L has the shape 

LP(X, i) = - E (Q.(X)P(X, t)) . (VI.3) 



For Langevin equations i.e. in (VI.2 1 Qi is replaced by Qi + rji where {rji} are Gaussian processes with < rjiit) >= 
and < rji{t)rij{t') >= 2AijS{t — t'), one has 

- - E + E A..^^(X,^) . (VI.4) 

i ah 

Let us introduce the eigenvalues {oik\ and the eigenfunctions {ipk\ of C: 

L^k afeV'fe • (VI. 5) 

Of course ipo — Pinv and ao = 0, and typically in mixing systems Reak < for k — 1, 2, .... Furthermore assuming 
that coefficient {gi, g2, ■■■} and {hi, /12, ...} exist such that functions ^(X) and /i(X) are uniquely expanded as 

g{X) = ^ gkM^) , H^) - E ^fcV^fc(X) , (VI.6) 

k=0 fe=0 

SO we have 

Cgj{t) = Y,9khk<ijl>e^^\ (VL7) 

k=l 

where Cg,f{t) —< g{'K{t))h{'K{t)) > — < .g(X) >< /i(X) >. For "generic" functions g and /, i.e. if they are not 
orthogonal to V'l so that 51 7^ and hi ^ 0, at large time the correlation Cg,f{t) approaches to zero as 

C„,(t)^e-/- , r..^. (VI.8) 

In some cases, e.g. very intermittent systems like the Lorenz model at r ~ 166.07, Reai = so the decay is not 
exponentially fast. 
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VII. APPENDIX C: LYAPUNOV EXPONENT IN DYNAMICAL SYSTEMS WITH NOISE 

In systems with noise, the simplest way to introduce the Lyapunov exponent is to treat the random term as a 
time-dependent term. Basically one considers the separation of two close trajectories with the same realization of 
noise. Only for sake of simplicity consider a one-dimensional Langevin equation 

dx dV{x) 

where r]{t) is a white noise and V{x) diverges for | a; |— > oo, like, e.g., the usual double well potential V = 

The Lyapunov exponent associated with the separation rate of two nearby trajectories with the same realization 
of r}{t), is defined as 

A,, = lim - \n\z{t)\ (VII.2) 

t— ^oo t 



where the evolution of the tangent vector is given by: 

dz d'^V{x{t)) 



it). (VII.3) 



dt 5a;2 

The quantity Ao- obtained in the previous way, although well defined, i.e. the Oseledec theorem (Bohr et al., 1998) 
holds, it is not always a useful characterization of complexity. 

Since the system is ergodic with invariant probability distribution P{x) — Cie^^'-^-'/'^^, where Ci is a normalization 
constant and C2 = one has: 

A, = limt^^ i In \z{t)\ = - limt^^ \ J* dl,V{x{t'))dt' 

(VII.4) 

= -Ci / dl^V{x)e-^^^y^^ ^x = -% /(a^V^(a;))2e-^(^)/^^ da; < . 

This has a rather intuitive meaning: the trajectory x(f,) spends most of the time in one of the "valleys" where 
—d'^^V{x) < and only short intervals on the "hills" where —d'^^V{x) > 0, so that the distance between two 

trajectories evolving with the same noise realization decreases on average. The previous result for the ID Langevin 
equation can easily be generalized to any dimension for gradient systems if the noise is small enough (Loreto et al., 
1996). 

A negative value of A^^ implies a fully predictable process only if the realization of the noise is known. In the 
case of two initially close trajectories evolving under two different noise realizations, after a certain time Ta, the two 
trajectories can be very distant, because they can be in two different valleys. For cr ^ 0, due to the Kramers formula 
(Gardiner, 1990), one has ^ e^^/*^^, where AV is the difference between the values of V on the top of the hill and 
at the bottom of the valley. 

Let us now discuss the main difficulties in defining the notion of 'complexity' of an evolution law with a random 
perturbation, discussing a simple case. Consider the ID map 

x{t +1) = f[x{t),t] + aw{t), (VII.5) 
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where t is an integer and w(t) is an uncorrelated random process, e.g. w are independent random variables uniformly 



distributed in [—1/2, 1/2]. For the largest LE A^, as defined in (VII. 2|, now one has to study the equation 

z{t + l)^ f[x{t),t]z{t), (VII.6) 

where /' = df/dx. 

Following the approach in (Paladin et al., 1995) let x{t) be the trajectory starting at a;(0) and x'{t) be the trajectory 
starting from a;'(0) — x{0)+6x{0). Let Sq = \Sx{0)\ and indicate by ti the minimum time such that |a:;'(Ti)— a:(Ti)| > A. 
Then, we put x'{ti) — x{ti) + dx{0) and define T2 as the time such that \x'{ti + T2) — x{ti + T2)| > A for the first 
time, and so on. In this way the Lyapunov exponent can be defined as 

X=-\n(^\ (VII.7) 



being t = '^Ti/M where A/" is the number of the intervals in the sequence. If the above procedure is applied by 



considering the same noise realization for both trajectories, A in (VII. 2 1 coincides with A^ (if Xa > 0). Differently, by 



considering two different realizations of the noise for the two trajectories, we have a new quantity 

K„=^ In , (VII.. 



which naturally arises in the framework of information theory and algorithmic complexity theory: note that K^/\n2 
is the number of bits per unit time one has to specify in order to transmit the sequence with a precision (5o, The 
generalization of the above treatment to A^-dimensional maps or to ordinary differential equations is straightforward. 



If the fluctuations of the effective Lyapunov exponent ^{t) (in the case of (VII. 5 1 ^{t) is nothing but In \ f'{x{t))\) 
are very small (i.e. weak intermittency) one has K„ = A + 0(ct/A) . 

The interesting situation happens for strong intermittency when there are alternations of positive and negative 7 
during long time intervals: this induces a dramatic change for the value of K^;. Numerical results on intermittent 
maps (Paladin et al., 1995) show that the same system can be regarded either as regular (i.e. A^r < 0), when the 
same noise realization is considered for two nearby trajectories, or as chaotic (i.e. K^j > 0), when two different noise 
realizations are considered. We can say that a negative Ao- for some value of a in not an indication that "noise induces 
order" ; a correct conclusion is that noise can induce synchronization. 
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Fig. 1: Lorenz-96 model: autocorrelation Cjj{t) (full line) and self-response Rjj{t) (+) of the fast variable Ukjit) {k = 3,j = 3). 
The statistical error bars on Rjj{t) are of the same size as the graphic symbols used in the plot. 
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Fig. 2: Lorcnz-96 model: autocorrelation Ckk{t) (full line) and self-response Rkk{t), with statistical error bars, of the slow 
variable Xkit) (fc = 3). 




Fig. 3: Lorenz-96 model: autocorrelation Cz^it) (dashed line) of the cumulative variable Zk(t) compared to the autocorrelation 
Ckkit) oi Xk{t) (full line). 
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Fig. 4: Loronz-96 model: time signal sample of the slow variable Xk{t) (k = 3) for the deterministic model (full line) and for 
the stochastic model (dashed line). For clarity, the two signals have been shifted from each other along the vertical a:xis. 
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Fig. 5: Loronz-96 model: PDFs of the cumulative variable Zk {k — 3), see definition in the text for the two cases, for the 
deterministic model (full line) and the stochastic model (dashed line). 
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Fig. 6: Loronz-96 model: PDFs of the slow variable Xk {k = 3) for the deterministic model (full line) and the stochastic model 
(dashed line). 



25 




0.5 1 1.5 2 2.5 

t 



Fig. 7: Lorcriz-96 model: autocorrelation Ckk{t) (full line) and self-response Rkk{t), with statistical error bars, of the slow 
variable Xkit) for the stochastic model. 
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Fig. 8: Lorenz-96 model: quadratic cross-response function Rfjf(t) for the deterministic model (full line), for the stochastic 
model when the slow variables evolve with the same noise realization for all components except one (dashed line), and when 
the slow variables evolve with a different noise realization for every component (dotted line). 
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Fig. 9: DW model with e = 1: time signal sample of the slow variable qo{t)- The ratio between fast and slow characteristic 
times is e ~ 0.1 (see text). 




Fig. 10: DW model with e = 1: time signal sample of the fast variable qi{t). 




Fig. 11: DW model with e = 1: autocorrelation Cii(t) (full line) and self-response Rii{t), with statistical error bars, for the 
fast variable gi. 
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Fig. 12: DW model with e = 1: Autocorrelation Coo(t) (full line) and self-response Roo{t), with statistical error bars, for the 
slow variable go- 
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Fig. 13: DW model with ? = 0.01. implying e ~ 10"'': autocorrelation Coo(i) (dashed line), self-response R(m{t), with 
statistical error bars, and the correlation function C{t) predicted by the FRR (full line) which is actually undistinguishable 
from the response. 
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Fig. 14: Comparison of the PDFs of the transition time te between the two climatic states for the DW model (full line) and 
the WNDW model (dashed line), for ?= 1 (e ~ 0.1). 




Fig. 15: WNDW model: time signal sample of the slow variable qo{t)- 
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Fig. 16: PDFs of the slow variable qo for the DW model with e = 1, i.e. e ~ 0.1 (full line), the WNDW model (dashed line) 
and the DW model with ?= 10~^, i.e. e ~ 10"'^ (dotted line). In the limit e — > 0, the PDFs of the deterministic model and of 
the stochastic model collapse. 




Fig. 17: WNDW model: autocorrelation C\]i,){t) (dashed line), self-response R^o{t), with statistical error bars, and the 
correlation function C(t) predicted by the FRR (full line). 
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Fig. 18: Quadratic cross-response function Kfj!(t) for the DW model (full line) and the WNDW model (dashed line). The 
growth rates of _R* (t) for the DW model are compatible with the two characteristic times of the system, while for the WNDW 
model 7?^j {t) quickly saturates in a very short time. 



